#MM <- 1000
pos.vec <- c(100,200,300,400,500,600,700,800,900,1000)
bias.mat.base<-MAE.mat.base <-MAE.mat <-RMSE.mat.base <-RMSE.mat <-bias.mat <- st.mat.base <- st.mat <- matrix(NA,length(pos.vec),7)
MAE.mean.base <-  ARMSE.mat.base <-MAE.mean <-ARMSE.mat <-avg.bias.mat <- matrix(NA,length(pos.vec),1)
xlabtxt <- "Number of \n Hand-Coded Documents"
k <- 1                      
#route <- "/nfs/home/D/dhopkins/condor/"
route <- "/nfs/projects/p/poliblog/poliblog.cvs/condor/output/"

for(i in 1:length(pos.vec)){
  MM <- pos.vec[i]
  txt1 <- paste(route,"resmatN",MM,".dat",sep="")
  txt2 <- paste(route,"actmatN",MM,".dat",sep="")
  txt3 <- paste(route,"tmatN",MM,".dat",sep="")

  txt4 <-paste("resmat",MM,"<- resmatH <-read.table(file=txt1)",sep="")
  txt5 <-paste("actmat",MM,"<-actmatH <- read.table(file=txt2)",sep="")
  txt6 <-paste("tmat",MM,"<-tmatH <- read.table(file=txt3)",sep="")

  eval(parse(text=txt4))
  eval(parse(text=txt5))
  eval(parse(text=txt6))

  txt1 <- paste(route,"resmatO",MM,".dat",sep="")
  txt2 <- paste(route,"actmatO",MM,".dat",sep="")
  txt3 <- paste(route,"tmatO",MM,".dat",sep="")

  txt4 <-paste("resmat",MM,"<- resmatJ <-read.table(file=txt1)",sep="")
  txt5 <-paste("actmat",MM,"<-actmatJ <- read.table(file=txt2)",sep="")
  txt6 <-paste("tmat",MM,"<-tmatJ <- read.table(file=txt3)",sep="")
  
  eval(parse(text=txt4))
  eval(parse(text=txt5))
  eval(parse(text=txt6))

  txt1 <- paste(route,"resmatP",MM,".dat",sep="")
  txt2 <- paste(route,"actmatP",MM,".dat",sep="")
  txt3 <- paste(route,"tmatP",MM,".dat",sep="")

  txt4 <-paste("resmat",MM,"<- resmatK <-read.table(file=txt1)",sep="")
  txt5 <-paste("actmat",MM,"<-actmatK <- read.table(file=txt2)",sep="")
  txt6 <-paste("tmat",MM,"<-tmatK <- read.table(file=txt3)",sep="")
  
  eval(parse(text=txt4))
  eval(parse(text=txt5))
  eval(parse(text=txt6))
  
  resmatI <- rbind(resmatJ,resmatH,resmatK)
  actmatI <- rbind(actmatJ,actmatH,actmatK)
  tmatI <- rbind(tmatJ,tmatH,tmatK)  
#  resmatI <- rbind(resmatH)
#  actmatI <- rbind(actmatH)
#  tmatI <- rbind(tmatH)  

  
  ###QI 1: bias for each category
  ####estimate - truth
  bias.mat[k,] <- apply(resmatI-tmatI,2,mean)
  bias.mat.base[k,] <- apply(actmatI-tmatI,2,mean)
  
  ###QI 2: sum of absolute bias
  #avg.bias.mat[k,1] <-mean(apply(resmatI-tmatI,2,mean))

  ###QI 3: standard deviation of each category
  st.mat[k,] <-apply(resmatI,2,sd)
  st.mat.base[k,] <- apply(actmatI,2,sd)
  
  ###QI 4: RMSE
  RMSE.mat[k,] <-sqrt(apply(resmatI-tmatI,2,mean)^2 + apply(resmatI,2,sd)^2)
  RMSE.mat.base[k,] <-  sqrt(apply(actmatI-tmatI,2,mean)^2 + apply(actmatI,2,sd)^2)
  
  ###QI 5: RMSE on average over categories
  ARMSE.mat[k,1] <- sqrt(mean(apply(resmatI-tmatI,2,mean)^2 + apply(resmatI,2,sd)^2))
 ARMSE.mat.base[k,1] <- sqrt(mean(apply(actmatI-tmatI,2,mean)^2 + apply(actmatI,2,sd)^2))
  
  ###QI 6: MAE, mean average error
  MAE.mat[k,] <- apply(abs(resmatI-tmatI),2,mean)
  MAE.mat.base[k,] <- apply(abs(actmatI-tmatI),2,mean)
  
  ###QI 7: MAE, sum of mean average error
  MAE.mean[k,] <- mean(apply(abs(resmatI-tmatI),2,mean))
  MAE.mean.base[k,] <- mean(apply(abs(actmatI-tmatI),2,mean))
  
  k <- k+1
}


###QI 1: Bias by Category

#postscript(file="/nfs/fs1/projects/poliblog/poliblog.cvs/paper/bias.eps",height=10,width=20)

xlabtxt <- "Number of Hand-Coded Documents"

pdf(file="/nfs/fs1/projects/dhopkins/readme/wds/figs/biasA.pdf",height=5,width=5)
par(pty="s")
#par(mfcol=c(1,2))
labs <- c("Extremely Negative","Negative","Neutral","Positive","Extremely Positive","Not a Blog","No Opinion")
lbs2 <- c("-2","-1","0","1","2","NB","NO")

for(i in 1:7){
  text1 <- "Nonparametric Estimator"
 plot(pos.vec,bias.mat[,i],main=text1,xlab=xlabtxt,ylab="Bias",type="l",cex=1.2,cex.lab=1.2,cex.main=1.6,cex.axis=1.2,ylim=c(-.1,.1),lty=i,xlim=c(0,1000),axes=F)
axis(2,pos.vec,at=c(-.10,-.05,0,.05,.10),labels=c("-.10","-.05",".00",".05",".10"),cex.axis=1.2)
axis(1,cex.axis=1.2)
par(xpd=T)
for(j in -4:4){
  lines(y=c(j/100,j/100),x=c(-40,-45),lty=1)
}
par(xpd=F)     
  abline(h=0)
  text(lbs2[i],y=bias.mat[1,i],x=28,cex=1.2)
  if(i %in% 1:6) par(new=T)
}
dev.off()

pdf(file="/nfs/fs1/projects/dhopkins/readme/wds/figs/biasB.pdf",height=5,width=5)
par(pty="s")
labs <- c("Extremely Negative","Negative","Neutral","Positive","Extremely Positive","Not a Blog","No Opinion")
lbs2 <- c("-2","-1","0","1","2","NB","NO")
for(i in 1:7){
text1 <- "Sampling Estimator"
  plot(pos.vec,bias.mat.base[,i],main=text1,xlab=xlabtxt,ylab="Bias",type="l",cex=1.2,cex.lab=1.2,cex.main=1.6,,cex.axis=1.2,ylim=c(-.1,.1),lty=i,xlim=c(0,1000),axes=F)
axis(2,pos.vec,at=c(-.10,-.05,0,.05,.10),labels=c("-.10","-.05",".00",".05",".10"),cex.axis=1.2)
axis(1,cex.axis=1.2)
par(xpd=T)
for(j in -4:4){
  lines(y=c(j/100,j/100),x=c(-40,-45),lty=1)
}
par(xpd=F)   
abline(h=0)
  text(lbs2[i],y=bias.mat.base[1,i],x=28,cex=1.2)
  par(new=T)
}
dev.off()


##QI 4


pdf(file="/nfs/fs1/projects/dhopkins/readme/wds/figs/rmsebycat.pdf",height=8,width=8)
par(pty="s")
                                        #postscript(file="/nfs/fs1/projects/poliblog/poliblog.cvs/paper/rmsebycat.eps")
mx2<-max(RMSE.mat)
plot(pos.vec,ARMSE.mat,main="",ylab="",type="l",xlab=xlabtxt,cex=1.6,cex.lab=1.8,cex.main=2.5,cex.axis=1.5,ylim=c(0,mx2))
par(new=T)
plot(pos.vec,ARMSE.mat.base,main="",ylab="Avg. Root Mean Squared Error",type="l",xlab=xlabtxt,cex=1.6,cex.lab=1.8,cex.axis=1.5,cex.main=2.5,lty=2,col="red",ylim=c(0,mx2))


dev.off()


